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1. Introduction 

Since tlie 19tli century potential theory has been growing continuously 
as gravitational and electromagnetic interactions and forces were derived 
from potentials satisfying Laplace, Poisson and Helmholtz equations. Some 
boundary value problems such as the Dirichlet problem and the Neumann 
problem, the electrostatic distribution of charges on conductors or the Robin 
problem can all be defined in terms of potential theory. For domains with 
sufficiently smooth boundaries the above problems use specific types of po- 
tentials like the volume potential, the single- and the double-layer potentials, 
the logarithmic potential for two-dimensional domains, etc. [l|, |2|. 

The Dirichlet and Neumann problems defined on domains with sufficiently 
smooth boundaries (i.e., a regular piecewise Lyapunov surface [3|) can be re- 
cast in integral equations which lead to compact operators on domain bound- 
ary: the Neumann- Poincare operator and its adjoint [l|, l2]- These methods 
are successfully applied in some practical and physical problems regarding di- 
electric heterogeneous systems [4'] like radio-frequency and microwave dielec- 
tric spectra of living cells |5| and plasmonic properties of metallic nanopar- 
ticles 0,0,9. 

Another problem which stems from potential theory is the equilibrium 
charge distribution on a conductor (the Robin problem) [9| and the implicit 
capacitance with applications in computational biophysics 10|, in scanning 



robe microscopy [ll|, ll2| , or in electrical charge storage in supercapacitors 



[l3|. The capacitance of an arbitrarily shaped body is directly related to 



hydrodynamic friction and to diffusion-controlled reaction rate [ij]. Thus a 



method of calculating the capacitance is based on "mimicking" the diffusion- 



controlled reaction rate on an object of arbitrary shape as random walks [15 
Another random walk method is defined just on the boundary and is given 
by the ergodic generation of the equilibrium charge distribution [16]. 

A standard procedure for solving the Laplace equation is the Finite Ele- 
ment Method (FEM) [l7| and the Boundary Integral Equation (BIE) method 
which leads to a finite element formulation as the Boundary Element Method 
(BEM) [l8|. In the FEM the entire domain is discretized by using elements 
and the associated basis functions, while in the BEM only the surfaces of the 
inclusions within the domain are used for discretization. The bounded and 
the unbounded domains are on equal footing in the BEM. On the other hand, 
the FEM solves the static and quasi-static electromagnetic field problems in 
unbounded domains by using one of the additional ingredients: infinite ele- 



ments 



19| . coordinate transformations 20|, or hybrid FEM/BEM methods 



2l| . Infinite domains can be elegantly and efficiently treated with the hy- 



brid Trefftz method, where, beside the usual elements, an additional set of 



functions is employed to treat singularities and infinities [22 



In this paper, by a BIE method |23l . |24| . |25| , we will define the capaci 



tance and the equilibrium charge of an arbitrarily shaped object in terms of 
the eigenvectors of the Neumann-Poincare operator and of its adjoint |l|, |2|. 
Compared to other methods the operator based BIE method provides directly 
the geometric dependence of capacitance that is incorporated in an eigenvec- 
tor norm of the Neumann-Poincare operator. This eigenvector corresponds 
to the largest eigenvalue. Accordingly, capacitance is obtained concurrently 
with other physical properties like localizedplasmon resonances of metal- 
lic nanoparticles that have the same shape p, |7|, |8|. Furthermore, operator 
perturbations on eigenvalues and on eigenvectors permit the perturbative 
estimation of capacitance and other physical properties. The reformulation 
of capacitance allows defining several scale-invariant shape factors that can 
also be readily used in the estimation of capacitance for arbitrary shapes. 
Applications regarding the capacitance of a general capacitor, the membrane 
capacitance of living cells, charge storage in supercapacitors, and the electric 
capacitance of molecular nanowires are discussed. 

The paper is organized as follows. The method is described in the second 
section. Then, the connection with a general capacitor is made in the fol- 
lowing section. Section 4 describes the numerical method, its accuracy and 
some applications. The conclusions are summarized in the last section. 

2. The operator approach in potential theory. Capacitance of a 
metalhc object 

We assume an arbitrarily shaped domain Q of volume V and bounded 
by the surface S. We denote the complementary set of Q hj Q'^ = ^^\Cl, 
where Q is the closure of Q and ^^ is the Euclidian 3-dimensional space. The 
following operators can be defined on E: 

1 f M (y) n (x)_(x - y) 

y| 



M[u] = - I ^-^^ ,J ^;3 ' dJ:{y), (1) 



x,yeS 

its adjoint. 



1 f V (y) n (y) ■ ( x - y) 

|x-y| 



^f*M=j^ / '•"::•"■.:: '' ''^(y). (2) 



x,yeS 

which is the Neumann- Poincare operator [2|, and the symmetric and non- 
negative Coulomb operator 

S[u] = j- I ^^rfS(y). (3) 

47r J |x-y| 
x,yes 

In Eqs. ([I])-(l2]) n is the normal vector to S. The operator M can be made 
symmetric via the Plemelj's symmetrization principle M^ S = S M }^. Thus, 
the operator M is symmetric with respect to the inner product defined by 
the symmetric and non-negative operator S 

{v\u), = (v\S[u]), (4) 



where ( | ) defines the standard inner product on L^ (S) and ( | )^ is the inner 
product determined by S. The Hilbert space L^ (S) is the vector space of 
all square- integrable functions defined on S. The standard inner product on 
L? (S) of two functions ui (x) and U2 (x) is defined as 

(mi I M2) = / ul {x) U2 (x) (is (x), (5) 

xGS 

where the star sign * signifies the complex conjugate of a complex number. 
The operators M and M"^ have the same spectrum bounded by the interval 
[-1/2 ,1/2] and the eigenfunctions Ui of M are related to the eigenfunctions 
Vi of M^ hj Vi = S [ui], which makes them bi-orthogonal, i.e. if M [ui] = Xi^i 
and M^ [vj] = Xj'^ji then {vj \ Ui) = 5ij [g, |2|, |23|]. In physical terms M [u] 
is related to the normal component of the electric field that is generated 
by surface charge density u. On the other hand, Vi is the electric potential 
generated on surface S by the charge distribution Ui. The largest eigenvalue 
of M and M^ is 1/2 irrespective of the domain shape and its corresponding 
eigenfunction vi of M^ is a constant function, i. e., Vi =constant on S, which 
comes from the following relation [2] 

constant f n( y) ■ (x - y) ^^ ^_^ _ 1 

|x-y| 



M' M = ^== ""' '^'" dS (y) = ± ., (6) 

47r ./ x-v 2 



yes 



Equation (Q can be interpreted in terms of the solid angle under which S can 
be seen from an arbitrary point located also on S. As we will see below, the 
companion eigenfunction ui of M is proportional to the equilibrium charge 
distribution on a conductor of shape determined by S. We note that the 
spectrum of M and ikf^ is invariant under the transformation x — )■ tx, t > 
(the scale invariance). In contrast, the spectrum of S is proportional to the 
linear size of Q. This can be easily checked in the case of a sphere where the 
eigenvalues of S are proportional to the radius of the sphere. Thus, because 
of S, Ui and Vi are not scale invariant. Therefore, the norms of Ui and Vi are 
proportional to the square-root and to the inverse of the square-root of the 
linear size of Q, respectively. The size-dependence of Ui and Vi is also reflected 
in the the bi-orthogonality conditions defined in Ref. |6|. Operators dl])-® 
can be used in the resolution of many physical problems. Dielectric spectra 
of living cells in microwave and radio- frequency regimes |5|, |23| or plasmonic 
properties of metallic nanoparticles [g, |7|, |8|, |2J, |25[ can be treated with the 



help of M and M^ whose inversion is carried out via the eigenfunctions Mj 
and Vi and the corresponding eigenvalues Xi- However, the above problems 
do not require normalized eigenfunctions because Ui and Vi come in pairs 
in the desired solutions p, la, |7|, |8|, |23|, |2J, |25[. Thus, there is no need for 



explicit calculation of S. Operator 5" is used whenever an explicit expression 
for the near-field properties like the near-field enhancement produced by 



localized plasmon resonances are needed [25|. Another example where 5* 
would be needed is the calculation of capacitance and of equilibrium charge 
on a metallic object. The Robin problem of finding the equilibrium charge 
distribution on a conductor of arbitrary shape can be cast into the operator 
form as |9| 

M [ur] = -ur (x) , (7) 

with the constraint / ur (x) rfS (x) = 1. The constraint (normalization 

xes 
condition) can be put in the following form 

(1 \ur) = 1 , (8) 

where 1 e L^ (S) is the constant distribution that has the value 1 on S. 
Equation ([71) is found from the jump formula obeyed by the derivative of 
the single-layer potential across E J2| and it has the obvious solution ur cx 
Ml. In other words, the solution of Robin problem is proportional to the 



eigenfunction of the largest eigenvalue of M . The constant value Vr of the 
electric potential generated by u^ is formally given by 



S [ua] = Vnl 



(9) 



and it is called the Robin constant, while its inverse is the capacitance C 
of the body bounded by E [l[. The left-side hand of Eq. ^ needs to be 
divided by dielectric permittivity e of the embedding medium if any physical 
situation is considered. Thus, from (|8]) and ([9]), in the general case of an 
outer medium with dielectric permittivity £, the capacitance is: 



C 



The quantity Wnj^ = \Ur 



1 


e 


'^« {un 


S [ur] ^ 


SM)/(2e) 


is just t 



(10) 



is just the electrostatic energy of 

the charged metallic body with charge density ur. The electrostatic energy 
Wuj^ of the equilibrium charge density ur is the energy minimum attained 
for the set of arbitrary surface charge densities that satisfy the constraint 
provided by Eq. (jH]). This is the Thomson theorem J26| and it is used as the 
definition for the equilibrium charge distribution in abstract mathematical 
terms [l|. The proof is as follows. For any surface density u one has the 

oo 

expansion u = Y2 o^iUi-, where Oj are the expansion coefficients of the surface 

density u. In addition, u obeys (1 |m) = 1, which fixes the first expansion 
coefficient since (1 \ui) = for i ^ 1 (the biorthogonality condition). One 
can prove that the constant ai is the proportionality factor between vi and 
the constant distribution 1, i. e., 1 = vi/ai. Thus 



u 



S[u]) = E ((^>i S[ajUj] 

JO oo 

Y^ a*aj {ui \vj) = Yl O'lO'j^ii 



oo 



El 1^ ^ I I 
\ai\ > \ai\ 



Ur 



i=l 



S[ur\ 



(11) 



which demonstrates the Thomson theorem. Moreover, the explicit expression 
of the equilibrium distribution is 



Ur = aiUi 



(12) 




Figure 1: (Color online) Schematic representation of a capacitor with the configurations 
of fields inside Si and outside S2. The dotted lines delimitate a Gaussian surface that is 
used further in the main text by invoking the Gauss law. 



and ai is related to the norm of f 1 by 

(fi |fi) = ||fi|| = al (1 |1) = a^A, (13) 

where A is the area of S. Finally, considering Eqs. (IT2l) and ( 1T3l) the capac- 
itance expression (ITOl) takes a simple and compact form 



C 



sA 

l^ill' 

|2 



(14) 



We have pointed out earlier that ||fi||^ is proportional to the linear size of 
O, therefore the capacitance itself is also proportional to the linear size of 
the body. Equation flT4|) shows explicitly both the geometric dependence of 
capacitance of an arbitrarily shaped object and the scale invariance of the 
shape factor C/\/47rA. The shape factor varies slowly with the conductor 
shape [27[ thus, one can say that ||fi|| is a slowly varying function of the 
conductor shape. In some cases like that of a cube |28i] the capacitance is 
difficult to compute, but as it is suggested by Eq. ([H]) operator perturbations 



can be used to estimate the capacitance of any object [29 



3. Capacitors and their capacitance 

The capacitance of an arbitrary object is defined and calculated with 
respect to infinity where the electric potential is supposed to vanish. In 



general, a capacitor consists of two separated conducting bodies. We consider 
a capacitor that is made of two closed and smooth surfaces Ei and E2 in which 
S2 encloses Si (Fig. [1]). The capacitor capacitance is defined as the total 
charge that is held on Si when the electrical potential is 1 on Si and on 
S2. This definition can be cast in a form in which the total charge on Si is 
given by the Gauss law or by the total electrostatic energy enclosed within 



the space between Si and S2 [26[. The electrostatic problem of a capacitor is 
to find a function that is 1 on Si and on S2 and obeys the Laplace equation 
in the space Qu between Si and S2: 

Am (x) = 0; X G Qu 

«(x) = 1; xe Si (15) 

M (x) = 0; X G S2. 

It is easy to see that, with the boundary condition taken from (llSp . the 
solution of the Laplace equation inside of Si and outside of S2 is the constant 
1 and the constant 0, respectively. Inside fli2 we seek a solution for ( ITSjl in 
the form of two single-layer potentials 

Ait J |x — y| 47r J |x — y| 

yeSi yGE2 

where fii and /i2 are the induced charge densities on Si and S2. Similar to 
([1]) we define on Si and S2 four operators Mn, M12, M21, and M22 as 

1 f f^j (y) n (x) • (x - y) 



M^AH] = J- / "^^-^^ ^ ' ,3 ^rfS(y), (17) 

47r J |x — y| 

xGE^yeSj 

with i,j = 1, 2. The equations obeyed by /ii and /i2 are 



(18) 



Mil [/^l] + M12 [fi2] = i/ii 
M21 [fil] + M22 [H = -|/i2, 

which say that the normal fields on Si from inside and on S2 from outside 
are zero. The solution of flTSj) is the solution of flTSj) up to a multiplicative 
constant. The multiplicative constant is fixed by the equations that set the 
boundary conditions of f lT^ : 

-^ll [^1] + ^12 [^2] = 1 /-.QN 

^2lN+^22[H=0, ^ ^ 

8 



with Sij being similar to (j2]) 

S.Au] = ^ I 0^\dJ:iy). (20) 

The capacitance of the capacitor is the total charge represented by charge 
density /ii and depends not only on operators like ([T]) and (^ but also on 
inter-surface operators flTTl) . In the special case when S2 is an equipotential 
surface determined by the equilibrium charge distributed on Si a compact 
formula based on operators ([I])- ([2]) can be deduced for capacitance. It is 
not hard to see that Eqs. ( ITSl) have as solutions the charge densities /xi and 
/i2, which are proportional to the equilibrium charge densities induced on Ei 
and S2, respectively. In addition, /ii and /i2 must have opposite signs. To 
determine fii and fi2, and the capacitor capacitance one needs also Eqs. (TT9|) . 
Thus, by integrating the first equation of ( !T9|) on Si and the second equation 
on S2 one obtains the following relations 

where Vi is the electric potential induced by fii on Si, V2 is the electric 
potential induced by fi2 inside S2 as well as on Si, and V12 is the electric 
potential induced on S2 by ^i. On the other hand, the total charge is 

Qi= J /ii(x)t/Si(x) = CiV^i (22) 



xGS 



on Si and 



Q2= I /i2(x)dS2(x)=C2y2 (23) 

XGE2 



on S2. We mention that Eqs. (1221) and ( 1231) are valid only whenever S2 
is one of the equipotential surfaces determined by an equilibrium charge 
distributed on Si. Equation f fT4|) provides the expressions of Ci and C2 that 
are the capacitances of Si and S2, respectively. Moreover 

Qi + Q2 = (24) 



which is the Gauss law that is obtained from the second equation of (ITS]) . 
Combining Eqs. ( 1?!]) . (12^ . (125]) . and (121|) we obtain the capacitor capacity 



C. 



cond 



(25) 



1 1 

Ci 6*2 

Particular examples in which Eq. (!25|) is directly applicable is the capac- 
itance of concentric spheres and of a coaxial cable. For a capacitor made 
of two concentric spheres the capacitance is Csph_cond = ^^3^? where Ri 
and i?2 are the radii of the two spheres with R2 > Ri- Having in mind the 
capacitance of a sphere as Csph = 4:7ieR, it is easy to check that Csph_cond 
has the form given by Eq. f l25|) . Moreover, the capacitance of a capacitor 
made of two confocal spheroids obeys also f l25|l . with Ci and C2 as having 
analytic expressions that are given in the standard textbooks of classical 
electrodynamics |26| . 




aspect ratio 



Figure 2: (Color online) Calculated capacitance per unit area versus aspect ratio for 
spheroids (black squares) and for cylidrical rods with different end-cap geometries: oblate 
spheroid (red triangles), sphere (green stars), and prolate spheroid (blue diamonds). It is 
assumed that the largest cross-sectional radius is 1 m. 



4. A Numerical method. Discussion 

As we have previously discussed, capacitance can be determined in numer- 
ical simulators used to calculate other physical properties like the localized 



10 



plasinon resonances in metallic nanoparticles p, |7|, |8|. The method is an 
operator based BIE method that calculates numerically the eigenvalues Xk 
and the eigenvectors Uk and Vk of M and M"^ , respectively. In order to have 



normed Uk and Vk one needs also to calculate the matrix elements of 5 [25 . 



The method belongs to the class of the spectral methods, where the functions 



of the basis set are defined globally rather than locally like in the FEM 30|. 
Details of the method for axisymmetric objects are given in Refs. 23|, |24 |. 
We have performed calculations for oblate and prolate spheroids as well as for 
cylindrical rods with different end-cap geometries: half of an oblate spheroid 
with 1/2 aspect ratio, half of a sphere, and half of a prolate spheroid with an 
aspect ratio of 2. To set the input data we have considered that the largest 
cross-sectional radius is 1 m for all the above geometries. The aspect ratio 
of any considered object is the ratio between the largest axial length and the 
largest cross-sectional diameter. The results are given in Fig. [21 where the 
capacitance per unit area is plotted. For the same aspect ratio the rods have 
larger capacitances than the spheroids. Moreover, for aspect ratios greater 
than 5 the capacitance of the rods do not depend any longer on the end-cap 
geometry. 

Our numerical calculations show a very good agreement with the ana- 



lytical results for spheroidal shapes which can be found in Refs. [26|]. The 
relative error is at least 5 x 10^^ with a relative small overhead of 25 func- 
tions in the basis and 96 quadrature points. We have also checked these 
results against the more sophisticated Trefftz based FEM calculations pro- 
vided by the multi-physics program ANSYS [31l]. The relative error in this 
case is around 10~^ (about 1% and two order of magnitude greater than the 
BIE results). In ANSYS the capacitance results obtained with the Trefftz 
method are much better than the results obtained by using infinite elements. 
On the other hand, we compared the BIE calculations with ANSYS Trefftz 
calculations for the rods and the two calculations are apart by less than 1%. 
An example of rods with finite length is that of capped carbon nanotubes 
which are used as electrochemical double-layer capacitors (supercapacitors) 



for energy storage [32j. The implementation of BIE for axisymmetric shapes 
has also shown to provide very accurate results of the depolarization factors 
which are related to the other eigenvalues of M and M^ [33] . 

A direct application of these calculations is the estimation of membrane 
capacitance in living cells. The shelled ellipsoidal model is one of the most 
common models of living cells in modeling dielectric spectroscopy experi- 



ments [3J]. In the ellipsoidal shelled model the shell designates the cell 

11 



membrane that is bounded by two confocal ellipsoids and it is practically 
non-conductive in radiofrequency. Living cells accumulates positive/negative 
charge on the outer/inner surface of the membrane, giving rise to a resting 
potential across the membrane. In dielectric spectra the charge on the mem- 
brane is responsible for the a-relaxation (usually below the frequency of 10 
KHz), whilst the dielectric mismatch between the cell and the surrounding 
medium gives rise to /9-relaxation in the MHz range of the radiofrequency 



spectrum J35j. The cell resting potential may be an input parameter in theo- 



retical analysis of both a- and /3-relaxations, thus the membrane capacitance 



needs to be estimated [36|, l37| . In the spherical model of living cells Prodan 
et al. 37| used the expression of a parallel-plane capacitor for the normalized 
capacitance (capacitance divided by area). It turns out that the formula of a 
parallel-plane capacitor is very close to that of a spherical living cell. The cell 
membrane is very thin with respect to the overall size of the cell, hence, for 
a spherical model of living cells, the membrane capacitance is approximated 
by 

Csyhjmemhr'^eA2ld = C2l{b-\), (26) 

where d = R2 — Ri, with Ri and R2 as the inner and outer membrane radius 
of areas Ai and A2, respectively. C2 is capacitance of the outer sphere and 
6 = R2/R1 > 1 obeys the condition (5 — 1) <^ 1. In the general case of 
arbitrarily shaped cells the parallel-plane capacitor formula cannot be used 
anymore although intuition may suggest otherwise. An analytical formula 
for the membrane capacitance of spheroidal livingcells is provided by the 
combination of capacitance formula for spheroids [26] and Eq. |25l 

In Fig. |3] there are presented two shape factors that are related to capac- 
itance and are scale-invariant. These shape factors can be straightforwardly 
utilized in approximate capacitance calculations for metallic object of vari- 
ous shapes. The first scale-invariant shape factor (= 74^'^/(27r^'^||fi|p)) is 
related to the area of the object (Fig. |3^) and it shows a relative shape 
insensitivity for aspect ratios less than 5. This shape factor has been used in 
isoperimetric inequalities to estimates the capacitance of objects with shapes 



close to the spherical shape [27[. In fact. Chow and Yovanovich noticed that 
this shape factor varies little over a wide range of shapes. 

The second scale-invariant shape factor related to capacitance is V^^^^/| |'i^i| | 
It is defined by using the volume of the object and shows shape insensitivity 
for long structures (Fig. [3]d). Thus for an aspect ratio greater than 5 the 
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Figure 3: (Color on-line) Scale-invariant shape factors versus aspect ratio. The shape 
factors are defined (a) by the area of the object and (b) by its volume. The spheroids 
are denoted by black squares and the cylidrical rods with different end-cap geometries are 
depicted by red triangles (oblate spheroid end-cap), by green stars (spherical end-cap), 
and by blue diamonds (prolate spheroid end-cap). 
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rods and the spheroids have almost the same volume related shape factors, 
which also vary little not only with the shape but also with respect to the 
aspect ratio. We note here that quantum capacitance of molecular nanowires 



is scaling with V^^^ rather than with the length of the nanowire 38| , thus for 
long structures the volume rather than the area of the object plays a greater 
role in determining the capacitance of an object. 

5. Summary and Conclusions 

We have presented a very compact formula of the capacitance for an 
arbitrarily shaped metallic object. Capacitance can be calculated with the 
help of the Neumann- Poincare operator, which together with its adjoint plays 
a central role in the resolution of elliptic partial differential equations like the 
Laplace equation. The formula of capacitance is stated simply as follows. 
The capacitance is direct proportional to the dielectric permittivity of the 
embedding medium and to the area of the object, and inverse proportional 
to the squared norm of the eigenfunction of the Neumann-Poincare operator 
with the largest eigenvalue. The operator approach on capacitance permits 
an elegant proof of the Thomson theorem which says that the electrostatic 
energy accumulated on a metallic body reaches its minimum when the charge 
distribution on the body is in equilibrium. In addition, the capacitance and 
the charge distribution on a metallic object is a byproduct obtained in BIE 
numerical simulators of localized plasmon resonances. 

We have shown also how the capacitance of a capacitor can be related 
to the individual capacitance of each surface of the capacitor. Thus, when 
the outer surface of the capacitor is an equipotential surface generated by 
the charging of the inner surface then the capacitor behaves like a series 
capacitor with the total capacitance as being the capacitance of the inner 
surface in series with the opposite (negative) capacitance of the outer surface 
of the capacitor. These results have been used in the analysis of membrane 
capacitance in living cells, where the parallel-plane capacitor model works for 
spherical cells but is not appropriate for membrane capacitance of arbitrarily 
shaped living cells. 

The capacitance formula allows us to define scale-invariant shape factors 
that can be used in the approximate calculation of capacitance. We have 
analyzed two scale-invariant shape factors. One of the shape factors employs 
the volume of the object and is more suitable for long shapes like rods or 
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wires. Alternatively, the other shape factor, which is defined in terms of the 
object area, is appropriate for objects with shapes close to a sphere. 

6. Acknowledgments 

This work was supported by a grant of the Romanian National Authority 
for Scientific Research, CNCS - UEFISCDI, project number PNII-ID-PCCE- 
2011 -2-0069. 

References 

[1] O. D. Kellog, Foundations of Potential Theory, Springer- Verlag, Berlin- 
Heidelberg-New York, 1967. 

[2] D. Khavinson, M. Putinar, H. S. Shapiro, Arch. Ration. Mech. Anal. 
185 (2007) 143. 

[3] N. Giinter, La theorie du potentiel et ses applications aux problemes 
fondamentaux de la physique mathematique, Gauthier-Villars, Paris, 
1934. 

[4] D. J. Bergman, Phys Rep. 43 (1978) 377. 

[5] D. Vrinceanu, E. Gheorghiu, Bioelectrochem. Bioenerg. 40 (1996) 167. 

[6] F. Ouyang, M. Isaacson, Philos. Mag. B 60 (1989) 481. 

[7] I. D. Mayergoyz, D. R. Fredkin, Z. Zhang, Phys. Rev. B 72 (2005) 
155412. 

[8] T. J. Davis, K. C. Vernon, D. E. Gomez, Phys. Rev. B 79 (2009) 155423. 

[9] J. G. Van Bladel, Electromagnetic Fields, Wiley-Interscience, Hoboken, 
New Jersey, 2007. 

[10] T. Simonson, Rep. Prog. Phys. 66 (2003) 737. 

[11] W. A. Hofer, A. S. Foster, A. L. Shluger, Rev. Mod. Phys. 75 (2003) 
1287. 

[12] S. Hudlet, M. S. Jean, C. Guthmann, J. Berger, Eur. Phys. J. B 2 (1998) 
5. 

15 



[13] P. Simon, Y. Gogotsi, Nature Mat. 7 (2008) 845. 

[14] J. B. Hubbard, J. F. Douglas, Phys. Rev. E 47 (1993) R2983. 

[15] J. F. Douglas, H. X. Zhou, J. B. Hubbard, Phys. Rev. E 49 (1994) 5319. 

[16] M. Mascagnia, N. A. Simonov, J. Comput. Phys. 195 (2004) 465. 

[17] C. Johnson, Numerical Solutions of Partial Differential Equations by 
Finite Element Method, Cambridge Univ. Press, Cambridge, 1987. 

[18] D. Poljak, C. A. Brebbia, Boundary Element Methods for Electrical 
Engineers, WIT, Boston, 2005. 

[19] P. Bettess, Int. J. Numer. Methods Eng. 11 (1977) 53. 

[20] D. A. Lowther, E. M. Freeman, B. Forghani, IEEE Trans. Magn. 25 
(1989) 2810. 

[21] S. J. Salon, J. D'Angelo, IEEE Trans. Magn. 24 (1988) 80. 

[22] Q. H. Qin, Appl. Mech. Reviews 58 (2005) 316. 

[23] T. Sandu, D. Vrinceanu, E. Gheorghiu, Phys. Rev. E 81 (2010) 021913. 

[24] T. Sandu, D. Vrinceanu, E. Gheorghiu, Plasmonics 6 (2011) 407. 

[25] T. Sandu, Plasmonics xx (2013) xxx, doi:10.1007/sll468-012-9403-z. 

[26] L. D. Landau, E. M. Lifshitz, Electrodynamics of Continuous Media, 
Pergamon, Oxford- New York, 1984. 

[27] Y. L. Chow, M. M. Yovanovich, J. Appl. Phys. 53 (1982) 8470. 

[28] C. O. Hwang, M. Mascagni, J. Appl. Phys. 95 (2004) 3798. 

[29] D. Grieser, H. Uecker, S. A. Biehs, O. Huth, F. Riiting, M. Holthaus, 
Phys. Rev. B 80 (2009) 245405. 

[30] J. P. Boyd, Chebyshev and Fourier Spectral Methods, Dover, New York, 
2001. 

[31] www.ansys.com 



16 



[32] J. Huang, B. G. Sumpter, V. Meunier, G. Yushin, C. Portet, Y. Gogotsi, 
J. Mater. Res. 25 (2010) 1525. 

[33] T. Sandu, J. Nanopart. Res. 14 (2012) 905. 

[34] K. Asami, T. Hanai, N. Koizumi, Jpn. J. Appl. Phys. 19 (1980) 359. 

[35] R. Stoy, K. Foster, H. Schwan, Phys. Med. Biol. 27 (1982) 501. 

[36] C. Prodan, E. Prodan, J. Phys. D: AppL Phys. 32 (1999) 335. 

[37] E. Prodan, C. Prodan, , J. H. Miller, Biophys. J. 95 (2008) 4174. 

[38] J. C. Ellenbogen, C. A. Picconatto, J. S. Burnim, Phys. Rev. A 75 (2007) 
042102. 



17 



